CERN-TH/99-356 
TTP99-55 
frep-ph/0001062 
December 1999 



ON 



Oh: 



S3 



Quarkonium momentum distributions 
in photoproduction and B decay 



o 
o 
o 

M. Beneke and G.A. Schuler 

P-i- 

<^ ■ Theory Division, CERN, CH-1211 Geneva 23, Switzerland 



and 



^ i S. Wolf 

\q . Institut fur Theoretische Teilchenphysik, Universitat Karlsruhe, 

D-76128 Karlsruhe, Germany 

O 
O 

Abstract 

Oh 



According to our present understanding many J/ifi production processes 
proceed through a coloured cc state followed by the emission of soft par- 
ticles in the quarkonium rest frame. The kinematic effect of soft particle 
emission is usually a higher-order effect in the non-relativistic expansion, 
but becomes important near the kinematic endpoint of quarkonium en- 
ergy (momentum) distributions. In an intermediate region a systematic 
resummation of the non-relativistic expansion leads to the introduction of 
so-called 'shape functions'. In this paper we provide an implementation 
of the kinematic effect of soft gluon emission which is consistent with the 
non-relativistic shape function formalism in the region where it is appli- 
cable and which models the extreme endpoint region. We then apply the 
model to photoproduction of J/ip and J /ip production in B meson decay. 
A satisfactory description of B decay data is obtained. For inelastic char- 
monium photoproduction we conclude that a sensible comparison of theory 
with data requires a transverse momentum cut larger than the currently 
used 1 GeV. 



PACS Nos.: 13.85.Ni, 14.40.Gx 



I. INTRODUCTION 



Inclusive charmonium production processes can be expressed in a factorized form, com- 
bining a short-distance expansion with the use of a non-relativistic QCD Lagrangian 
(NRQCD) [I]. The short- distance expansion works best for total production cross sec- 
tions, provided the expansion parameter v 2 (of order of the typical velocity squared of 
the quarks in the bound state) is small enough. It follows that contrary to prior belief 
many charmonium production processes such as production in hadron-hadron collisions 
at large transverse momentum || and at fixed target 0, and in B meson decay [§-0], 
are actually dominated by production of a coloured cc state, followed by a long-distance 
transition to charmonium and light hadrons ||. 

The theoretical prediction of charmonium energy distributions is more delicate. 
A long-standing problem for the NRQCD factorization approach concerns the z- 
distribution in inelastic J /if) photoproduction, where z = Ej/^/E^ is the quarkonium 
energy fraction in the proton rest frame. The colour octet contributions to this quantity 
grow rapidly near z = 1 |9|||, in conflict with observation [10|, unless the NRQCD matrix 
elements that normalize the colour octet contribution are made rather small.Q 

One of the physical origins of this discrepancy is as follows: the fragmentation of 
the coloured cc state into J /if) occurs via the emission of gluons with small momentum 
fractions of order v 2 . Because the momentum of these gluons is small compared to 
the momenta involved in the hard subprocess that creates the cc state, it is neglected 
in leading order in the short-distance expansion (in v 2 ); the fragmentation into J /if} 
is described by a single number (the 'NRQCD matrix element'). This is adequate for 
total production cross sections, but it is not for distributions in the kinematic region, 
where the charmonium carries nearly maximal energy. In this region, the J/ if) energy 
distribution is evidently sensitive to the energy distribution of the soft emitted gluons. 
In particular, we expect that the J /if) energy distribution should fall to zero, rather 
than grow, near the point of maximal energy, if the J/ if) is produced via a colour octet 
state, since the emission of gluons with momentum much smaller than their typical one 
is rather unlikely. 

The inadequacy of a leading-order treatment of the short-distance expansion, and 
the necessity to account for the kinematics of soft gluon emission, is even more evident 
for J I ip production in B meson decay. The leading order partonic short- distance process 



* There may be other difficulties for the NRQCD factorization approach, which we do not 
discuss in this paper. For a long time, transverse polarization of J /if) produced in hadron- 
hadron collisions at large transverse momentum |11-13| has been regarded as the crucial test of 
the theoretical framework. If recent indications from CDF of no polarization are confirmed 



by higher statistics data, this may indicate a problem with factorization, as suggested in [15|, 
or it may imply large spin-symmetry violating corrections. 
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b — ► (cc)g results in cc pairs with fixed (maximal) energy, in stark contrast to the broad 



energy distribution observed ||16|| . The broad energy distribution of multi-body final 
states has to be attributed to soft gluon emission and to the Fermi motion of the b quark 
in the B meson. 

Technically speaking, the velocity expansion of the short- distance process breaks 
down near the kinematic endpoint of maximal charmonium energy |T^,|T^], because 



higher-order terms in the small parameter v 2 are compensated by inverse powers of 
small kinematic invariants. Such a breakdown of the short- distance expansion is not 
specific to quarkonium production in the NRQCD approach, but occurs quite generally 
for inclusive processes, for example in deep-inelastic scattering as Bjorken x — > 1 or in 
semi-leptonic or radiative B decays jnj. When the quarkonium carries a fraction (1 — e) 
of its maximal energy, where e is small, the inclusiveness of the process is restricted by 
the small phase-space left for the emission of further particles. The process is then also 
sensitive to the fact that the physical phase space is limited by hadron kinematics, while 
the calculation of short- distance coefficients is carried out in terms of partons. The 
short- distance expansion reacts to this non-inclusiveness by exhibiting terms of order 
{y 2 /e) k . In some cases one can sum the leading terms in v 2 /e to all orders and express 
the quarkonium production cross section as a convolution of a non-perturbative 'shape 
function' with a partonic cross section. The shape function leads to a smearing of the 
energy spectrum. The shape function formalism is analogous to a leading twist approx- 
imation, and is appropriate for e ~ v 2 . In this intermediate region the framework of 
the NRQCD factorization approach is still valid, reorganized by a partial resummation 
of the velocity expansion. However, in the extreme endpoint region, e <C the twist 
expansion also breaks down. 

The leading twist expressions for several energy distributions have been derived in 
fl8|| . But since the shape function is non-perturbative and essentially unknown, no quan- 
titative analysis has been performed. It is the aim of this paper to explore the kinematic 
effect of soft gluons in the fragmentation of a coloured cc pair quantitatively. In particu- 
lar, we will be interested in the question whether folding the short-distance cross section 
with a shape function can indeed account for the observed ^-distribution in J/ip photo- 
production. The emission of soft gluons with energy of order m c v 2 in the quarkonium 
rest frame cannot be computed perturbatively and we have to model it. Our ansatz for 
the soft gluon radiation function will be guided by simplicity. The important feature of 
the model is that it incorporates the kinematics of soft gluon radiation together with 
reasonable assumptions on the typical energy scales involved. The ansatz bears some 
similarities with Fermi motion smearing [ED] and, in particular, the ACCMM model [2I| 
for semileptonic B decays. Since the precise form of the energy distribution near the 
endpoint depends on the ansatz for the shape function, our results do not constitute the- 
oretical predictions. However, as we shall see, a satisfactory description of B decay data 
can indeed be obtained with a reasonable ansatz for the shape function. A further cross 
check is provided by applying the same shape function to the J/ ip energy distribution in 
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photoproduction. This however, turns out to be more problematic. 

The paper is divided as follows: Sect. 2 is 'theoretical'. We define the model and 
derive the equation that describes the convolution of the short-distance process with 
the shape function for a general production process. We also show that the model 
is equivalent to a specific form of the NRQCD shape function in the region where a 
leading twist approximation is valid. To illustrate the formalism, we consider the limit 
m c v 2 ^> Aqcd, in which charmonium is a Coulomb bound state. We rederive NRQCD 
factorization for this specific case and compute the shape function in this limit. 

The ansatz for the non-perturbative shape function depends on a few model parame- 
ters. In Sect. 3 we apply the model to the J/ip momentum distribution in B — > J/ip + X 
and tune the parameters of the model to the observed momentum distribution. In Sect. 4 
the more complicated (and more interesting) case of inelastic photoproduction is consid- 
ered. 



II. SHAPE FUNCTION MODEL 

In this section we derive the general expression for the smeared quarkonium energy 
distributions on which the applications to B decay and photoproduction will be based. 
To motivate our approach and to make more explicit contact with the formalism of [H,0 , 



we begin by considering the production amplitude for quarkonium in the Coulomb limit, 
and with emission of a single soft gluon, before generalizing the expressions to the case 
of interest. In the last subsection we return to the Coulomb limit and compute the shape 
function in this limit. This provides us with an idea of the form of the shape function 
in a controlled, although unrealistic limit. 



A. Factorization and the shape function in the Coulomb limit 

Inclusive charmonium production proceeds in two stages [0: first a pair of nearly on-shell 
and co-moving charm quarks is created in a hard process in which typical momenta are 
of order 2m c (or larger, if there is another hard scale) in the charmonium rest frame. The 
nearly on-shell cc state then fragments into charmonium via emission of soft particles 
with energy and momentum of order m c v 2 in the charmonium rest framed Schematically, 
the differential cross section is expressed in the factorized form 

(27r) 3 2p fl -^ = flux • / dPStpi,^-] (27r) i 5(p R + J2k J + J2^- p in) 
d Pr J 4 i 



^The energy scale for these particles is set by the small velocity v that characterizes the non- 
relativistic charmonium bound state and the typical virtuality (m c v) 2 of the nearly on-shell c 
and c quark. See also the discussion below. 
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Pn 

FIG. 1. Diagrammatic representation of the amplitude that leads to Eq. ( |2,1| ). 

■ H(P in , P, h, l 2 , Pi ) S(p R , P, h, l 2 , kj), (2.1) 

where dPS\pi,kj] denotes the phase space measure for the sets of hard (p^ and soft 
(kj) particle momenta and H and S refer to the hard and soft parts of the amplitude 
squared, respectively. See Fig. [I] for a graphical representation and further explanation 
of not at ion. 

To define the hard and soft parts in ( |2.1| ) accurately, we use the amplitude for the 
process — > J/tpgg, relevant to inelastic photoproduction, as an example. It is also 
instructive to take the limit m c v 2 3> Aqcd, where v is now of order a s (m c v) and Aqcd 
is the strong interaction scale. We call this the Coulomb limit, because the charmo- 
nium bound state is perturbatively calculable in this limit and the dominant binding is 
through the Coulomb force. The Coulomb limit is much stronger than the non-relativistic 
limit. While charmonium and bottomonium are non-relativistic (v 2 <C 1), they are not 
Coulombic (m c v 2 ~ Aqcd) in reality. In particular, the NRQCD matrix elements, which 
are usually taken as non-perturbative parameters, can be perturbatively calculated in 
the Coulomb limit, up to corrections suppressed by powers of Aqcd / (m c v 2 ) . 

A particular contribution to the jg — > J/ipgg amplitude is shown in Fig. |2|. The 
corresponding squared amplitude is the sum of terms where both gluons are hard or both 
gluons are soft or one of them is hard and the other is soft. The hard-soft term is the 
most interesting one for inelastic photoproduction through the colour octet mechanism 
and we focus on it first. The other two terms will be briefly discussed later. 

Suppose the gluon with momentum px in Fig. ^|is hard and the gluon with momentum 
k is soft. On-shell soft gluons in NRQCD can have energy of order m c v and m c v 2 p2 



(called 'soft' and 'ultrasoft', respectively, in p2fl ). However, gluons with energy of order 



*In an abuse of notation, in the figure H and S refer to the hard and soft part of the amplitude, 
rather than the amplitude squared. The nearly on-shell heavy quark propagators that connect 
the hard and soft part in the figure should be considered as part of S. See below. 
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Px 



FIG. 2. A contribution to the 75 — > J/ipgg amplitude discussed in the text. 

m c v cannot be radiated over the time-scale l/(m c v 2 ) and do not appear as final state 
particles in the scattering amplitude.^ Consequently, the scale of k is m c v 2 . This is 
important, because this will set the scale for the energy of soft gluon emission in our 
model parametrization later. 

In Fig. |2] we included (dashed lines) the instantaneous exchange of (Coulomb) gluons 
with energy of order m c v 2 and momentum of order m c v. If this exchange occurs between 
nearly on-shell heavy quark propagators with off-shellness of order m c v 2 , it is not sup- 
pressed by the small coupling constant, because the total contribution from each gluon 
is of order a s (m c v)/v ~ 1. However, if one of the heavy quark propagators is far off- 
shell, Coulomb exchange represents an ordinary higher-order correction to the amplitude. 
Hence we can neglect gluon exchange to the left of the gluon with momentum px- The 
gluon ladder 'between' the emission of the gluons with momentum px and k, respectively, 
cannot be neglected, but it is summed into the Coulomb Green function G c (p,p'] E), the 
Green function for the Schrodinger equation with the (leading order) Coulomb potential. 
The Green function is related to the quark- ant iquark scattering amplitude for a quark- 
antiquark pair with (small) relative three momentum 2p into a quark-antiquark pair 
with (small) relative three momentum 2p' with total non-relativistic energy E. Likewise 
the gluon ladder to the right of the gluon with momentum k is summed and contained 
in the bound state wave function. For a 3 S\ state, such as J/i/j, the bound state wave 
function in the bound state rest frame is given by 



§More technically, because the interaction with a gluon with energy of order m c v sends the 
heavy quark propagator off-shell, a subgraph with energy and momentum of order m c v in the 
amplitude squared has no cc + ng cut, as would be required for a non-zero contribution to the 
73 — ► J/ipgg amplitude. Rather such a subgraph can be expanded into a series of instantaneous 
interactions, which contribute to the potential between the heavy quarks. 
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*(p fl ,A;p) = ^2M~ R - 



Sab 



(2.2) 



^ v/2 



where 



ip(p) 



8y^ 7 5 / 2 

(p 2 + 7 2 ) 2 



(2.3) 



and 7 = m c Cpa s /2. M R is the quarkonium mass, 5 a b refers to colour (with N c = 3 
the number of colours and (7p = (iV 2 — l)/(2iV c )), cr* is a Pauli matrix and e*(A) the 
polarization vector of the quarkonium. 

With these remarks one of the two (symmetric) hard-soft contributions to the am- 
plitude can be written as 



denotes the hard sub-amplitude with the on-shell spinors for its external heavy quark 
lines with momentum P/2 + l\ and P/2 + I2 removed. We also introduced the vector P, 
defined as P = (2m c , 0) in the J/ip rest frame, the relative momentum q = (h — h)/2 
and E{p R + k) = p° R + k° — 2m c = —m c {CFCt s ) 2 /'i + k°. The binding energy at leading 
order has to be kept in the last expression, because it is of the same order as k°. For 
later use we define I = l\ + 12, the vector that describes the motion of the cc pair in the 
J/ip rest frame. Note the kinematic relation P + I = p R + k. 

The amplitude is not yet in a factorized form, because the hard sub-amplitude still de- 
pends on q and / and its spin and colour indices are entangled with those of the remaining 
part of the amplitude. As described in [Q, we can perform a spin and colour decom- 
position that disentangles the two parts of the amplitude. We then expand the hard 
sub-amplitude in the small momentum q, which amounts to an expansion in derivative 
operators and a decomposition in orbital angular momentum. As a matter of principle, 
we could also expand the hard sub-amplitude in /. However, since it is / that occurs 
in the phase space constraint and that is related to the terms in the short-distance ex- 
pansion, which we intend to sum to all orders, we do not perform this expansion. The 
spin and colour decomposition, and the expansion in relative momentum q, results in 
the following expansion of the amplitude squared: 




(2.4) 
ccg) 



1-4(73 -> J l^99)\ 2 = J2 Pr « \Ml9 -> ccg)] Pr' n \Ml9 -> ccg)* 



n 




- tr [T n (q) iG c (q, p + fe/2; E(p R + k)) V(k; p) V(p R , A; p)] 




- tr [T' n {q) iG c {q, p + fe/2; E(p R + k)) V(k; p) V(p R , A; p)f • (2.5) 
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Here T n (q) is a matrix in spin and colour indices and a polynomial in q. The operator 
Pr n is also a matrix in spinor and colour indices and extracts the appropriate Taylor 
coefficient of expansion of A(jg — > ccg) in q. The quantity Pr n A(jg — > ccg) is q- 
independent, but still depends on I. In conventional NRQCD terms, the sum over n 
corresponds to intermediate cc pairs in different angular momentum and colour states, 
and also to higher dimension operators in each intermediate channel. The previous 
equation can be written as the product of a hard and soft part, 

\A{ 19 - J/^gg)\ 2 = £ H n (P m ,P,l, Px )S n (p R ,P,k), (2.6) 

n 

where the soft part S n is given by the last two lines of Q2.5p . H n and S n are still coupled 
through the relation P + 1 = p R + k, so we introduce 1 = / dH S(p R + k — P — I). Adding 
the phase space integration over px and k, we recover the differential cross section in a 
form similar to (12.11): 



r dH r 

= E Jj^M flux J dFS M (2ir) 4 6(P + l + p x -P m )H n (P m ,PJ,p x ) 

■ J dPS[k] (2n) 4 S(p R + k-P-l) S n (p R , P, k), (2.7) 

where a(cc[n])(l) refers to the short-distance part and F n (l) to the soft part. The ex- 
pansion in local operators appropriate to integrated cross sections [|l| is recovered after 
expansion of H n (P in , P, l,p x ) in I- In leading order, we then identify 

J dPS[k}S n (p R ,P,k) (2.8) 

with the NRQCD matrix elements defined in 

Before continuing let us discuss as an example the angular momentum and colour 
projection for the case of an intermediate cc pair in a 3 Si, colour-singlet state, at lowest 
order in the expansion in q. In this case Pr simply sets q to zero in the hard sub- 
amplitude and T(q) carries no q-dependence. The correctly normalized spin and colour 
projection is 

r » ( " )s>r » w - 75k' ff,0 75k' ff ' (2 - 10) 

where the trace includes a colour trace and the projection of the hard amplitude is 
written in a covariant form. Let us check that (2^) together with the projection ( p.lOj ) 
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do indeed reproduce the colour singlet NRQCD matrix element. In leading order the 
transition 3 S[^ — > 3 Si does not require gluon emission. Hence 



dPS[k]S n (p R ,P,k) 



E 



d 3 p tr(<x^(p K ,A;p)) 



(27T) 



M, 



B 



2m r 



(2.11) 



where we used that in the leading order approximation Mr pa 2m c . 

Note that F n (l) in (|2.7|) defines a more general object than the shape function in 
which is a function of only one variable l-\ 



lo + l z or l . The definitions of [IS] would 



be reproduced, if we could neglect the other components of I in the short-distance part. 
We shall discuss later, after generalizing ( |2.7p to the emission of more than one gluon, 
under what conditions this is justified. 

Up to now we considered the contribution of the diagram in Fig. to J/ip pho- 
toproduction, when one of the two emitted gluons is hard and the other is soft. The 
contribution from two hard gluons is part of the next-to-leading order correction to the 
short- distance part of the colour-singlet intermediate state. The contribution from two 
soft gluons smears out the contribution from the diagram with no gluon emission, which 
is concentrated at z — 1 and zero transverse momentum. It also contributes to the 
endpoint of the energy spectrum, but can be eliminated with a transverse momentum 
cut sufficiently large compared to several hundred MeV. Experimental measurements of 
inelastic J/ip photoproduction usually imply such a cut. 



B. The general case 

We now extend the previous discussion in the following way. We consider a general, 
inclusive charmonium production process (cf. Fig. |I|) 

Initial state (P in ) cc[n] + X(pi) J/ip(p R ) + X(p { ) + Y(kj), (2.12) 

where the cc pair is in a certain colour and angular momentum state n, X denotes a col- 
lection of hard particles, and Y a collection of soft particles emitted in the fragmentation 
of the cc pair. 

Since m c v 2 ~ Aqcd, the coupling to soft gluons is large and the emission of multiple 
gluons is not suppressed. Hence the emission of soft gluons is better described as the 
emission of a soft colour multipole field, which carries away a total momentum k = J2j hj 
and which has the correct quantum numbers to effect the transition from cc[n] to J/ip. 
Hence we define 

$ n (k;p R ,P) = J dPS[k 3 }(27r) 4 6(k-J2k,)S n (p R ,P,k :j ), (2.13) 
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where S n (p R , P, kj) is the generalization of the soft sub-amplitude that appears in (|2.7| ) 
to the emission of more than one soft gluon. With this definition the generalization of 
( p.7|) is given by 



(2vr) 3 2^ 



da 



E 



dH 



a(cc[n)){l) ■ F n {l) 



d*p R t^(2vr) 4 
J2 f-Sii flux / dPS[ Pi ] (2n) A 5(P + l + J2Pi~ p in) H n (P in , P, l,Pi 
dk 2 d 3 k 

" ' (27r) A 6(p R + k-P-l)$ n (k;p R ,P) 



2tt (2vr) 3 2A; 



(2.14) 



As above the differential cross section is factored into a short- distance and a soft part. 
In higher orders in the strong coupling, this would require careful subtractions to define 
both parts properly. We will be working only with cases, where the lowest order, tree 
approximation to the short- distance part is assumed. Then the factorization is trivial, 
as in the example of the previous subsection. 



There is an additional assumption implicit in writing ( 2.14 ), which concerns the 
validity of NRQCD factorization in general |l|, not only its generalization to spectra. 
The assumption is that the transition from the cc[n] state to J / ip occurs via emission of 
gluons rather than by absorption from the surrounding 'partonic medium'. Of course, if 
n is a colour octet state the emitted gluons must interact with the remnant process to 
form colour neutral hadrons; the NRQCD approach assumes that the process of colour 
neutralization is suppressed by powers of Aqcd/^c and can be formally ignored, if we 
consider v and Aqcd/^c as independent parameters such that AQ CD /m c <«<Cl. On 
the other hand, absorption would violate factorization explicitly, since its details depend 
on the environment created by the specific production process. Despite the fact that this 
issue affects most quarkonium production processes, it has rarely been addressed in the 
literature, with the exception of fT5(| . We will not dwell on this issue further and take 
factorization for granted. (The empirical fact that the NRQCD matrix elements are 
approximately universal, including hadronic collisions, may support this assumption.) 
However, an investigation of this point would certainly be useful. 



1 . Derivation of the smeared spectrum 



We now bring ( |2.14D into a more useful form. We make one additional simplification, 



which is adequate to the two applications which we consider in this paper. The simpli- 
fication is that there is only a single, massless hard particle in the final state. Then the 
set of momenta Pi consists of only p x , and p\ = 0. 

It is often convenient to refer explicitly to the quarkonium rest frame defined by 
p R = P = rather then the centre-of-mass frame defined by Pi n = 0. In the following 
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non-invariant quantities will refer to the quarkonium rest frame. For example, in {j)r — 
P) - l = {Mr — 2m c ) Iq, lo refers to the zero-component of I in the quarkonium rest frame. 
We define the z-direction as the direction of P — Pj„ in the quarkonium rest frame and 
in the centre-of-mass frame. With this unconventional definition of the z-direction in the 
centre-of-mass frame the boost from the centre-of-mass to the quarkonium rest frame is 
in the ^-direction and the transverse components defined with respect to this axis are 
invariant. 

We use the two 5-functions in ( 2.14j) to integrate over p x and k. Then define l± = 



lo ± l z and write 

dH dl + dlodlj 



(2vr) 4 32vr 4 
The ^-function left over from the second 5-function in (2.14) fixes 



(2.15) 



l 2 ± = (M R - 2m c f - 2{M R - 2m c )l + l + {2l - l + ) - k 2 . (2.16) 
The result of these manipulations is 

(2vr) 3 2^ ^ 



with 



d 3 p R 

r dk 2 dd> 1 

]T / —dl + dl 07 ^ S{A)-Mx-H n {P m ,P,l,p x ) ■ —<5> n {k; P R,P), (2.17) 

n J Z7T Z7T 47T 

A = {2m c - P m _) (2m c - P in+ + /+) + (2/ - 1+) {2m c - P in+ ) 

-{Mr - 2m c ) {Mr - 2m c - 2l ) + k 2 (2.18) 

and px = Pin — (P + I), k = P + I — pr. Furthermore, we have the constraints k > 0, 
p x ,o > 0, k 2 > and l]_ > 0. 

Any ansatz for the function Q n {k;pR, P) that we will be using will be independent 
of the azimuthal component <p oil. Hence we need only the azimuthally averaged short- 
distance part: 

H n {P m ,P,l,p x ) = J ^H n {P m ,P,l,p x ). (2.19) 

The remaining ^-function can be used to integrate over / + . Then we use k = 2m c — 
Mr + Iq as integration variable instead of lo and define 

a EE P m+ -Mr, (3 = P in . - Mr. (2.20) 

This leads to the final result 
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— J dk flux • ff n (P in , P, l,p x ) ■ _ $ n (k;p R , P). (2.21) 

" (a 2 +fc 2 )/(2a) 

Recall that a, /3 and ko are defined in the quarkonium rest frame. 

The integration limits are obtained as follows: inserting the constraint ( j2.18|) A = 
on l + provided by the last 5-function into l\ > with l\ given by (|2.16|) , we find the 
condition 

k 2 - a(2k - a)] [k 2 - (3(2k Q - (3)] < 0, (2.22) 

in addition to ko > and k < (a + /3)/2, which follows from px,o > 0. Now note that 
a(3 = (Pi n — pr) 2 > and that k > implies a + (3 > 0. Hence a and /3 are both 
positive. Now a — (3 = 2P in z . In the quarkonium rest frame the z-axis is defined by the 
direction of —Pin- This implies 

(3 > a > 0. (2.23) 



Eq. Q2.22j ) admits two solutions. The physical one yields the limits on the /co-integration 



in ( 2.21|) . The upper limit on the fc 2 -integral then follows. Note that k > and 



k Q < (a + (3)/2 are then respected automatically. 

Eq. fl2.21p is the main result of this section and we will use it later to obtain the J / ip 
energy spectra in B decay and photoproduction. Recall that flux • H n (P in , P,l,p x ) is 
just the ordinary, projected cc production cross section that enters familiar applications 
of NRQCD factorization with the only difference that the cc pair is produced with 
momentum P + I rather than P, and that an average over the azimuthal angle of / in 
the quarkonium rest frame is performed. This means that the invariant mass of the cc 
pair is given by M 2 5 = (P + l) 2 = (p R + k) 2 = M 2 R + 2M R k + k 2 >M 2 R rather than Am 2 c 
as in the conventional partonic calculation. This kinematic difference can make a large 
numerical effect. 

The radiation function & n (k;p R ,P) is defined by ( |2.13| ). Roughly speaking, it rep- 
resents the probability squared that a soft gluon cluster with energy k in the J/ip rest 
frame and invariant mass k 2 is emitted from the cc pair in the transition cc[n] — > J/i[>. 
We consider it as a non-perturbative function. We will make an ansatz and try to de- 
termine some of its parameters from existing data. In the Coulomb limit, the function 
&n(k;p R , P) could be computed as indicated previously. However, we shall not assume 
this limit for charmonium. 



2. The shape function limit 
As mentioned above, the function 
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F n {l) 



dk 2 d 3 k 
2tt (2vr) 3 2A; 



(2n) 4 6(p R + k-P-l)$ n (k;p R ,P) 



(2.24) 



defined in ( p. 14 ) is different from the shape function introduced in ||18|| . The shape 



functions introduced there correspond to a systematic resummation of enhanced higher 
order corrections in the NRQCD velocity expansion. Eq. ( |2.21| ) goes beyond such a 
systematic resummation. We now show that ( [2.14| ) and (|2.21|) are equivalent to the 
results of 0] in the region of applicability of the latter, up to non-enhanced higher 
order terms in the velocity expansion. 

We are concerned with energy spectra in a variable z. For quarkonium production in 
the decay of a heavier particle with mass m, we define z = 2Pj n • p R /P 2 n . The maximal 
value of z is z max = 1 + M^/m 2 , assuming that all other particles in the final state are 
massless. (In reality these will be pions; we neglect the small pion mass.) For quarkonium 
production in two-to-two collisions, a{pi) + b(p 2 ) — > J/ip + . . ., we define z = 2p 2 -p R /P 2 n . 
For example, in ^yp collisions p 2 is the momentum of the struck parton in the proton. 
The maximal value of z is z max = 1. 

Consider the ^-spectrum in the region z max — z of order w 2 <C 1, but z max — z not 



much smaller than v 2 . This is the region in which the shape function formalism of 
applies. We introduce px = J2i Pi i n ( |2.14j) and use the first 5-function to integrate over 
p x . This leaves a ^-function with argument 



[(P - P in )_ + L] [(P - P in )+ + 1+) - I] 



Px- 



(2.25) 



Using the definitions of z, it is easy to see that in the endpoint region p x and P — P in 
become nearly light-like. With our definition of the z-axis (P — Pj„)+ becomes small, of 
order m c v 2 (but not much smaller), while (P — Pi n )- remains of order m c .Q p\ has to 



be of order m 2 v 2 or smaller. All components of / scale as 



since Mr — 2m c and all 



components of k are of this order. It follows that the dependence of ( |2.25| ) on Z_ and 
l± can be dropped. Furthermore, the formalism of [18[] assumed that the dependence of 
the hard cross section H n (P in , P,l,p x ) on / can be neglected, since it is not related to 
enhanced higher order terms in the velocity expansion. As a consequence, we can pull 
the /_- and /^-integrations through to the second line of ( |2.14| ). The result then takes the 
form of a partonic differential production cross section convoluted with a shape function 
in l + , provided we identify the shape function defined in |18[ with 



F n (i + ) = ]T (o\^r nX m + Y)W + YW+ - ^+)(x t r»|o) 

Y 



dl_d 2 lj 
2(2vr) 4 



dk 2 d 3 k 
2tt (2vr) 3 2A; 



[2Ti)%p R + k-P-l) $ n (k;p R , P) 



(2.26) 



**A11 other large scales that the process may involve are treated as order m c . 
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This shows that ( |2.14| ) is consistent with the operator formalism of |I8| in the region of z 



where the operator formalism applies. Eqs. (|2.14j) and ( |2.21|) extrapolate this formalism 



into the extreme endpoint region z max - z C t) 2 . Since there is no correspondence with 
a systematic resummation of the velocity expansion in the extreme endpoint region, this 
extrapolation should be considered as a model. This is again analogous to energy spectra 
in semileptonic B decays []nj]. 

It is instructive to recover the consistency with the shape function formalism directly 
from ( 2.21|) . In the region z max — z ~ v 2 , we may approximate ( |2.18| ) by 



A « (2m c - P in _) (2m c - P in+ + 1+) = (a + M R - 2m c - l + ){(3 + M R - 2m c ). (2.27) 



This implies that the upper integration limits in Q2.21| ) are replaced by infinity.^ We can 



then re-introduce 1 = / dl+ S(l + — a — [Mr — 2m c ]) and factorize ( 2.21]) into a convolution 
over the hard cross section times the shape function ( 2.26|) . 



3. Form of $ n (k;p R ,P) 

Eq. ( |2.26| ) implies that the moments of the shape-function are related to the usual 
NRQCD matrix elements. For example, integration over l + results in 



/ l^Y Fn{l) = J2^I dk ^ k l- k2 ®n{k;p Rl P) = (0 J n ^), (2.28) 



oo oo 



(2tt) 4 ny > (2tt 







where (O^J^) is the conventional NRQCD matrix element for an intermediate cc pair in 
an angular momentum and colour state n. This could in principle be used to determine 
the overall normalization of & n (k;p R , P) from the known NRQCD matrix elements. 

In practice this is problematic. The phenomenological values of the NRQCD matrix 
elements are determined from integrated quantities in leading order in the velocity ex- 
pansion in a given channel n. On the other hand, if we compute the same integrated 
quantities from the spectra obtained with ( 2.21]) , they contain higher order terms in the 



velocity expansion, for example related to the fact that the invariant mass of the cc pair is 
always larger than the quarkonium mass M R . Since v 2 is not small, the integrated quan- 
tities can be quite different, if the normalization condition ( |2.28| ) is imposed. Another 
way of saying this is that the phenomenological values of the NRQCD matrix elements 
would be quite different from the commonly accepted ones, if the theoretical prediction 



TTThis is consistent with a ~ m c v 2 and (3 ~ m c in the shape function limit, such that 
af3 ~ m 2 v 2 and (3 + k 2 / [5 ~ m c , i.e. both upper limits are parametrically larger than the 
typical values of the integration variables k 2 ~ m 2 v 4 , and k$ ~ m c v 2 , respectively. 
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used to obtain them contained higher order terms in the velocity expansion. As a conse- 
quence we are forced to tune anew the overall normalization to the measured integrated 
spectra. We will return to this point below in the context of specific applications. 

The radiation function Q n (k;pR,P) is non-perturbative. Similar in spirit to the 
ACCMM model P]J for semileptonic B decays, we assume a simple functional ansatz 
for phenomenological studies: 

$ n (k;p R ,P) = a n ■ \k\ b " exp(-k 2 /A 2 n ) ■ k 2 ex P (-k 2 /A 2 n ). (2.29) 

The exponential cut-off reflects our expectation that the typical energy and invariant 
mass of the radiated system is of order A n ~ m c v 2 ~ several hundred MeV. Since the 
pattern of soft gluon radiation may depend on the cc state n, the parameters a n , b n 
and A n can differ for different states. The three parameters of the ansatz could be 
determined from the first three moments of the shape function. In practice this is not 
possible, not only because of the problem mentioned above, but also because the NRQCD 
matrix elements with derivatives to which the higher moments are related are not known 
phenomenologically. 

In later applications, we will need the radiation functions for the three colour octet 
states n = , 3 P (8) , 3 Sf } • We assume that 

b^] = 2, &[ 3 P (8) ] = b[ 3 S{ 8) ] = 0, (2.30) 

A[^ 8) ] = A[ 3 P (8) ] = A, A[ 3 sf } ] = cA. (2.31) 

The choice of ft^Sg ] = 2 is motivated by the fact that the gluon coupling for a Ml 
magnetic dipole transition from a ^Sq 8 ^ to J/ ip is proportional to the momentum of the 
gluon. Furthermore, the transition from 3 S{ 8) to Jji) occurs through two El electric 
dipole transition, which suggests that the average radiated energy and invariant mass is 
larger than for the single Ml and El transition in the other two cases. We fix c = 1.5; 
the effect of this somewhat arbitrary choice will be discussed in the context of specific 
applications. Of course, since soft gluon emission is non-perturbative for charmonium, 
the arguments that lead to these choices are at best indicative in any case. 

C. Computation of the shape function in the Coulomb limit 

In the following we compute the radiation function in the Coulomb limit m c v 2 ^> Aqcd, 
a s (m c v) ~ v for n = ^q 8 '* , 3 Pq 8 ^ to obtain an idea of the form of this function in a 
controlled limit. Since this limit is unrealistic for J/ip, the reader interested only in the 
application of the formalism presented above may jump directly to the next section.^ 



"The calculation is similar to a calculation reported in [23|. However, in this work the cc 
pair in state n is described by a Coulomb wave function just as J/tj). This substitution does 
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FIG. 3. Left: one of four NRQCD diagrams which contribute to the Ml (El) transition 
from an ^Sq ( S Pq ) state, specified by Tn\ to J/ip. Right: an example for a double El 
transition from a 

35(1,8) 

state. 

We begin with the chromo-magnetic dipole transition ccpjSg ] — > J/ip + g. With 
emission of one gluon ( J2.13[ ) simplifies to 

SpSfp^P) = 27r5(k 2 )S[ 1 S^](k;p R ,P). (2.32) 

Furthermore, S^Sq ](k; p R , P) is normalized to the conventional NRQCD matrix element 
according to (|2.8|), i.e. 

Ij02k° S ? S PK k '>P*> P ) = W^o))- ( 2 - 33 ) 

The non-relativistic quark-gluon vertices are classified according to their velocity 
suppression. The leading spin-flipping interaction is provided by the chromo-magnetic 
interaction vertex —g s /(2m c )(cr x k) with k the outgoing gluon momentum as in Fig. [3[ 
The diagram on the left hand side of Fig. gives (cf. ( |2.4j) ) 

glCp f d 3 k f k{kj\ f d 3 p d 3 p' d 3 q d 3 q' 
~8^f J (27r) 3 2fc° { ij ~ ~kF) J (2tt) 3 (2tt) 3 (2tt) 3 (2tt) 3 

• \ £ tr(e(A) ■*{* x k)i) tr(e*(A) ct(ct x (-fc));) i/,(p) ^{p') 

■ iG c (p + fc/2, q + fe/2; E(p R + k)) iG c {q' + k/2,pf + fc/2; £( P/? + fc)), (2.34) 

with as given by (p73|) and i?(p.R + fc) = P° R + k° - 2m c = -m c (C F a s ) 2 /4 + A; = 

— K 2 /m c . Eq. ( p.34|) can be simplified, because the gluon is ultrasoft with energy and 



not correspond to the NRQCD definition of a colour octet operator or the corresponding shape 
function, in which the cc pair is local and all intermediate states with the quantum numbers n 
are allowed, and described by the full Coulomb Green function. 
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momentum of order m c v 2 ~ m c a 2 , while p, p', q and q' are of order m c v ~ wi c o; s . 
Dropping small terms in the arguments of the Coulomb Green function (as we have 
already done when defining E(pn + k)), performing the traces and accounting for an 
identical contribution from the other three diagrams not shown in Fig. |3], we obtain 



S[ l Sn(k; PR ,P) = ^fk 2 imW, (2.35) 

m c 



ric(8)i 



2 



where 



/N 8) P) = J^0 i - 3 G c (q,p--K 2 /m c )^p). (2.36) 
To compute this integral, we switch to coordinate space, 

I[ 1 S$ l) ](k) = fd 3 xG c (x,0;-K 2 /m c )'ip(x), (2.37) 



use V^ 33 ) = y^/ne yx (7 = m c Cpa s /2), gained by Fourier transformation of ( p. 3D , and 
the following representation for the coordinate space Coulomb Green function 



G c (x, y; -K 2 /m c ) =J2( 2l + !) ^i{x, V\ -«7™c) Pi(x-y/(xy)), (2.38) 
z=o 

where x, y denote the modulus of a?, t/, ^(2) the Legendre polynomials and 



s=0 



Here L^ 2 ' +1 ^(2;) refers to the Laguerre polynomials and the parameter A is defined such 
that the Green function corresponds to the Green function in the potential 

^( r ) = _ A ^i. (2.40) 

r 

Hence A = 1, if the intermediate cc pair propagates in a colour singlet state, and A = 
— 1/ {2N c Cp) = —1/8, if it propagates in a colour-octet state, which is what we need here. 
Only the / = component of the Green function contributes to the integral ( [2.371) . The 



remaining radial integration over Laguerre polynomials is easily executed as an integral 
over the generating function 



e -zu/(l-u) l _ 



T = $> S L^) (2.41) 



==There is a misprint in the first reference of [24], which is corrected in Eq. (18) of the second 
reference. 
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with subsequent expansion in u. Then, summing over s, and introducing the dimension- 
less variable 

,1/2 

Z = fc/7 = 1 



Ak 



m c C'jpQ 2 S/ 



(2.42) 



the result is 

I[ l S {8 \k) 



( vr7 )l/2 ( z 2 



~ s(s-l/z) 

l) 2 ^i s-X/z Vl + z 



m. 



( 7T7 )l/2 ^2 _ X 
42 



1 + (A-1) 



z + 1 



_ (l _ 2 p x (-A/z, 1, 1 - A/z; (1 - z)/(l + *))) 



(2.43) 



with 2-^1(0, c; z) the hypergeometric function. Let us check the power counting: with 
7 ~ m c v and fc , k { ~ m c t> 2 , we obtain /[ 1 So ](&) ~ (m c /v) l l 2 and, from ( |2.33| ), (|2.35| ), 
(C^So)) ~ a s rr? c v' ' . This agrees with the velocity power counting of |]]. The additional 
a s arises, because we consider the weak coupling limit. 

The chromo-electric dipole transition cc[ 3 P ] — > J/ip + g is computed along similar 
lines. We have 



S[ 3 P (8) P;p R ,P) 



3mf. 



where 



d ~ 

— G c (x,y; -K 2 /m c ) 



and the normalization is given by 



d 3 k 



S[ 3 P (8) P;p R ,P) 



WP))> 



(2.44) 



(2.45) 



(2.46) 



(2vr) 3 2fc° " 

The derivatives in ( |2.45| ) come from the factor p in the electric dipole vertex —ig s {p + 
p')/(2m c ) ~ (—i)g s p/m c . In this case only the I = 1 component of the Green function 
survives the y — > limit and the angular integration. The result is 

4m c 7 3 / 2 z 3 ~ (s + l)(s + 2)(s + 3) (\ - z N 
37I 1 / 2 (z + lY±l 



iM 8) ](k) 



s=0 



s + 2-X/z 



l + z 



3/2 



(z + 1) 3 



2(1 + z)(2 + z) + (A - 1)(5 + 3z) + 2(A - if 



+ 



+ 



4z(l + z)(z 2 - A 2 ) 
(1-z) 2 

l + z)(z-A) 



iF 1 {-\/z, 1, 1 - A/z; (1 - z)/(l + z)) - 1 



(2.47) 
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FIG. 4. Dependence of A: /(4vr 2 ) S n (k;p R , P) for n 
gluon energy k$. The parameters are chosen as m, 



600 



lo(8) 



800 



Sq and n 



3d(8) 



on the emitted 



ml ~ a s mlv 7 , 



Velocity power counting gives (0${ 3 Po))/ 
the standard counting. 

The dependence of k /(4:ii 2 ) S n (k;pR,P) for n 



1.5 GeV, a s = 0.4, A = -1/8. 

which is again consistent with 



lo(8) 



Sq 0> and n = 3 Pq on the energy 



k of the emitted gluon is shown in Fig. |j. The input parameters are chosen as m c = 
1.5 GeV, ot s = 0.4; A = —1/8 for a colour octet matrix element. Both dependences are 
smooth and mainly reflect the asymptotic behaviours at small and large gluon energy. In 
particular the suppression of the curve at small k Q is a consequence of the structure 
of the magnetic dipole vertex. 

According to the normalization conditions ( |2.33| ) and ( 2.46| ) the integration of the 
two curves gives the value of the conventional NRQCD matrix elements. The result 
depends strongly (see the discussion below) on the cut-off on the integration range for 



k . Choosing the cut-off between 300 MeV and 600 MeV, we findQ 

(Ot^CSo)) = (0.07 - 0.61) • 10" 4 GeV 3 
(O J 8 ^( 3 P ))/m 2 c = (0.07 - 0.22) • 10~ 4 GeV 3 . 



(2.48) 



Although these numbers may be insignificant, because the assumption m c v 2 3> Aqcd 
necessary to obtain them, is not valid for charmonium, it is interesting to note that the 
matrix elements come out one to two orders of magnitude smaller than the phenomeno- 
logical values, determined from fitting colour-octet subprocesses to experimental data ||. 
This suggests either a large non-perturbative enhancement of the matrix elements - such 



***These numbers, in particular the one for 3 Pq, depend sensitively on a s . (0$^ (^Pq)) / m 2 
increases rapidly as a s increases. 
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as the presence of a gluon condensate to which the soft gluons can couple - or the pos- 
sibility that the phenomenological values of the matrix elements effectively parametrize 
other corrections to the production processes not related to soft gluon emission (such as 
higher order short-distance corrections). 

The behaviour of the soft function S n (k; pr, P) at large k deserves further discussion. 
First, we observe that the calculation by itself does not provide an intrinsic cut-off for 
large k$. This should not be expected, since at the level of perturbative radiation the 
ultraviolet behaviour of the soft function joins smoothly to the infrared behaviour of 
the short- distance part. A well known example of this occurs in P-wave production 
the logarithmic infrared behaviour of the coefficient function of (0* ( 3 Po)) matches the 
logarithmic ultraviolet divergence of (Og( 3 Si)). 

Inspection of S[Sq ](k;pji,P) shows that we obtain a quadratically ultraviolet di- 
vergent matrix element (Og^ ( 1 Sq)) , which seems to contradict the conventional wisdom 
that this matrix element is scale-independent at leading order. However, the conventional 
wisdom is derived from the use of dimensional regularization. If a hard cut-off on the 
gluon energy is used, the colour octet 1 Sq operator mixes into the colour singlet 3 5*i oper- 
ator through a quadratically divergent term.^ This corresponds to a infrared finite, but 
quadratically infrared sensitive contribution to the coefficient function of (O^ ( 3 5i)), 
consistent with the over-all v 4 suppression of (Og^^So)) relative to (0(^ ( 3 S\)). In di- 
mensional regularization, the quadratically infrared sensitive term is attributed entirely 
to the short-distance coefficient and the quadratic divergence in (Ot^fSo)) is set to 
zero. 

In case of S[ 3 PQ 8 ^](k;pn, P) we find a linear divergence for (Og^Po)). The inter- 
pretation of this divergence requires a more careful discussion of the fco-integral and the 
integrals over p and p'; it will not be presented here. 

In the ansatz (|2.29|) we have added a cut-off on ko by hand in the form of an expo- 
nential fall-off for ko ^> A n . We interpret this ansatz as a 'primordial distribution' for 
the radiation of non-perturbative gluons, which eventually is modified by perturbative 
evolution. This is similar to the assumption that intrinsic transverse momenta of the 
proton's constituents are bounded. Perturbative radiation violates this assumption and 
leads to the evolution of parton distributions. A similar ansatz is also implied by the 
ACCMM model or in shape functions for semileptonic B decays in general. 

Finally, we comment on the transition from a colour octet or a colour singlet 3 S\ state 
to J /if). This presents a more complicated case, since - besides the contribution with 
no gluon emission for the colour singlet state - the leading term requires the emission 
of two gluons, see the right hand side of Fig. [| In coordinate space this requires the 
evaluation of integrals of the form 



t't The dimensions work out correctly, because the two chromo-magnetic dipole vertices pro- 
vide two powers of l/m c . 
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iy[ 3 <Sf' 8) ] ~ Jd'xd'yG^-y^O^ipR + h)) 

x ^G c (x, y; E(p R + k 1 + fc 2 ))) (^0*0 J , (2-49) 

which we shall not pursue. If we were only interested in the limit of small loop momenta 
k = k\ + &2, we could expand the Green functions for k® <C 7 2 / m c first and integrate 
afterwards over x and y. We would then find the same small-fco behaviour as in the case 
of/[ 3 PoP). 



III. MOMENTUM SPECTRUM IN B -> J/V>X 

In this section we apply the formalism developed in the previous section to the J /if} 
momentum spectrum in the semi- inclusive decay B — > J/ipX. The leading partonic 
decay process is very simple, resulting in J/if> with fixed momentum, but the hadronic 
decay spectrum is modified by fragmentation of the cc pair, which is the main concern 
of this paper, and by bound state effects on the b quark in the B meson. Both will be 
taken into account in the following. 

We start by recapitulating the partonic result for b —>■ cc[n] + q. We then implement 
the fragmentation of the cc pair according to our shape function ansatz and obtain the 
J /if} momentum distribution in b quark decay. We regard this distribution as input 
distribution for the ACCMM model, which accounts in a simple but satisfactory way 
for the effect of Fermi motion of the b quark inside the B meson. The resulting J /if) 
distribution in B meson decay is then boosted to the CLEO frame and compared to 
CLEO data. The aim of this comparison is twofold: first we show that smearing of the 
spectrum due to fragmentation of the cc pair is essential to describe the CLEO data. 
Second we use these data to determine the shape function model parameter A. Assuming 
universality of the shape function over the whole kinematic domain, we will then turn to 



J/ip photoproduction in Sect. |V|. Results for the J /if) momentum distributions already 
exist in the literature, including colour octet production |2"5|J215|| . However, only Fermi 
motion effects are taken into account there. We will briefly compare our results with 



those of p5|,|26[] at the end of this section. 



A. Energy distribution in b quark decay 

The underlying partonic process of a B meson decay into J/if> and light hadrons is 
b — > cc[n] + q (q = {d, s}). Since the cc pair is treated as a single particle kinematically a 
leading order calculation of this process results in a fixed value for its energy (momentum) 
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rather than in a real spectrum. Defining^ z = 2E c5 / rrib as the energy fraction of the cc 
pair in the b quark rest frame, the "spectrum" is 

^k = T c - c 5(l + r ] -z), (3.1) 
dz 

where rj = Am 2 c /m1 for massless light hadrons in the final state. In a purely partonic 
calculation one may identify 2m c with the J/ip mass and m b with the B meson mass. 

At leading order in the non-relativistic expansion the cc pair has to be produced in 
a colour singlet 3 S\ state. This term coincides with the colour singlet model and has 
been computed long ago [p7| , p8[| . At relative order v 4 ~ 1/15 in the non-relativistic 
expansion, J / ip can also be produced through cc in 1 S^\ 3 Pj 8) , 3 S[ S) colour octet states. 
These formally subleading contributions are enhanced by a factor of about 15, by which 
the short- distance structure of the AB = 1 weak effective Hamiltonian favours the 
production of colour octet cc pairs in the b — > ccq transition. These additional terms can 
be comparable or even larger than the colour singlet term |^-[7[. They are the ones of 
interest in this paper, since the radiation of soft gluons in colour octet cc fragmentation 
has a large kinematic effect on the observed J/if> momentum spectrum. In comparison, 
fragmentation effects in the colour singlet channel are order v 4 suppressed relative to 
the total colour singlet rate and therefore negligible. Hard perturbative corrections to 
the colour singlet [0,^] and colour octet J7| production processes are also known. They 



enhance the colour octet channels moderately. Within the present limitations of the 
shape function ansatz we must neglect these perturbative corrections for consistency. 
The partonic production spectra for the cc[n] states of interest read 

= ^-P- H n (m b , 2m c ) 6{l + V -z), (3.2) 
dz 2m b 87T 

where 

nr<2 IT/ |2 rrj 4 

H n (m b , 2m c ) = C^f[n](r,) (3.3) 



and the process-specific functions f[n}(rf) are given by 

f[ 3 S[\r ] ) = (l- V )(l + 2r ] ), (3.4) 

f?S[ 8) ](v) = 1(1-^(1 + 277), (3.5) 



f[ 1 Sl ) 8) }(r ] ) = 9 -(l-r ] ), (3.6) 
f[ 3 PP]( V ) = 9(l-r ] )(l + 2 V ). (3.7) 



Win this section "hatted" quantities refer to the b quark rest frame. 
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Note that the colour octet matrix elements are not part of the hard amplitudes, but 
included in the normalization of the radiation function <& n (k), see ( |2.28|) . In case of the 
P wave contribution, the normalization refers to (O^Pq)) /m 2 and the corresponding 
factor 1/m 2 is also extracted from if [ 3 pj 8 ^](mb, 2m c ). As mentioned above we neglect 
QCD corrections and also small corrections due to penguin operators. The Wilson 
coefficients Cng] of the effective operators in the weak AB = 1 Hamiltonian are related 
to the usual C± by 

C w (/j) = 2C + 0j)-C_(Ai), 

q 8] (/i) = C + (/i) + C_(//). (3.8) 
At leading order, as appropriate to the present analysis, 

a s (M w 



C±(aO 



-7£ 0) /(2/3o) 



(3.9) 



with 7±' 1 = ± 2(3 =F 1) and /3q = 11 — 2nf/3. In ( |3.3|) the notation C[i )8 ] implies C[i], if 
n is a colour-singlet state, and C[g\, if n is a colour octet state. Note that CL/Cf-n ~ 15 
at /i ~ mi,. 

We now implement cc fragmentation for the colour octet production channels. Notice 
that the partonic amplitude squared has no azimuthal dependence, hence H n (nib, 2m c ) = 
H n (mb, 2m c ) in the notation of ( [2.19| ). Furthermore, we need the light cone components 
of the incoming momentum P in = (nib, 0) m the J/ip rest frame to get a and (3 of ( |2.20| ). 
We find 

a = JT^n - \p R \) -M R , (3= ^(E R + \p R \) - Mr. (3.10) 

M R Mr 

The index "R" now refers to J/ip. To complete the implementation we have to fix 
the ambiguity in treating the kinematic effects in the hard production amplitudes H n . 
Strictly speaking the shape function formalism allows us to ignore the dependence of the 
hard production process on the vector /, since it does not lead to singular contributions 
near the endpoint, if the hard matrix element is not singular at the endpoint. On the 
other hand, the invariant mass of the cc pair is kinematically given by 

MUk) = (p + I) 2 = (PR + k) 2 = M 2 R + 2M R k + k 2 , (3.11) 

where k is the four momentum of soft radiation in the J/i/j rest frame. We adopt the 
convention that 2m c in the partonic matrix element is replaced by M c5 (k) everywhere, i.e. 
even when it does not arise kinematically, but through internal charm quark propagators. 
This convention is consistent with the shape function formalism in the shape function 
limit, but is arbitrary otherwise. It has the advantage of incorporating the physically 
expected effect of reducing the short-distance amplitude, because of the need to create 
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a heavier cc pair as compared to a purely partonic picture. The only exception to the 
convention is the factor l/(2m c ) in ( |3.3| ), which comes from the normalization of the 
cc state. It should be replaced by l/M R . Eq. ([2.211) , specialized to the J /if) energy 
distribution in b quark decay, is then: 

a/3 2 03 2 +fc 2 )/(2/3) 

dTJj^ rdk i 1 ^ M« 

dE« 4tt» - 7 27r (a2+fe i /(2a) 2m 6 87rm 6 | P/? | 

with a, /3 from ( HTTP ). 



B. Normalization difficulty 



We assumed up to now that the radiation function $„(&) is normalized according to 
( p.28| ). This implies that as A of Q2.29D tends to zero the integral over dT/dE R of Q3.i2j ) 
equals the integrated partonic rate with m c = Mr/2. 

Consider now the integral T n (A) of the spectrum ( |3.12| ) with fragmentation (for a 
specific production channel n) at small A and expand in A. To make things simpler 
put k = in the hard matrix element H n . Then integrate over Er or, equivalently, 
z = 2En/mb, and perform a change of variables from z to a. Then note that for small A, 
one can set the upper limits of the k 2 and k integrations to infinity up to exponentially 
small corrections in A, given the ansatz ( p.29|) . Then exchange the k 2 and ko integration 
with the a integration to obtain 



f n (A) 



M, 



R 



16(2vr) 4 m fe 



H n (m b , M R ) / dk 



k + y/k 2 -k 2 

dk & n (k) I da 



Vk 2 



ko-^Jk 2 -k 2 



dz 



da 



(3.13) 



Now introduce the average 

1 



«/»» = 



oo oo 
2 



(2tt)3 ( W) J 



dk 2 / dk ^k 2 -k 2 ^ n (k)f(k), 



(3.14) 



defined such that ((l)) n = 1 according to the normalization condition ( 2.28|) . Eq. ( p,13| ) 
is then rewritten in the form 



f„(A) 



1 I — 77 



2m b 8ti 

where 77 is now defined as M^/m 2 and 



H n (m b ,M R ) (0^)-r n (A), 



(3.15) 



r„(A) 



M 



R 



2(1-77) 



U2 _ h.2 



k +^Jk 2 -k 2 

da 



k -^k 2 -k 2 



dz 



da 



)>»• 



(3.16) 
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Hence we obtain the partonic decay rate with m c = Mr/2 up to the factor r„(A). To 
evaluate r n (A) in an expansion in A we observe that 



dz 



CO 



w R [wrw - ") = ik I 1 - " + S ( " + 1] tjk) ) (3 ' 17) 



da 

can be expanded under the integral. The result is 

r » (A) = 1 + (- 2 « w R »■ + « »• + (si J ) • < 3 - i8 > 

The averages can be done using (|2.29|) with a n fixed by ( |2.28|) ; they scale with definite 
powers of A as follows from the form of ( |3.14j) . With 77 = 0.416, the result is 

/ * \ , f 4.761 A f 12.971 / A \ 2 . 1 . 

where the upper number refers to b n = in (|2.29| ) and the lower one to b n = 2. For 
A pa 300 MeV this implies large corrections to the integrated rate. Since A ~ m c v 2 , this 
must be interpreted as large higher order corrections in the velocity expansion, which 
are not taken into account in the usual leading order NRQCD analysis. This means 
that enforcing the normalization condition (|2.28|) underestimates the data, because the 
matrix elements on the right hand side of (|2.28|) have been obtained without these large 
higher order corrections. 

The effect is in fact even larger than indicated by ( |3.19| ), because we keep the k 
dependence of the hard matrix element and M c5 (k) is always larger than Mr. As an 
indication of this effect we can compute the average 

4mf 2 = {{M c - C {kf))n « Ml(l + j^j (3.20) 

which implies an effective charm quark mass of about 1.8 GeV rather than m c = 1.5 GeV 
which is usually adopted in partonic NRQCD calculations. 

When the implicit /c-dependence of the partonic matrix element H n is taken into ac- 
count, the numbers given in ( p.l9|) change. However, the observation that v 2 corrections 
are large is generic. 



C. Fermi motion effects 

We now convert the spectrum (|3.12) in b quark decay into a spectrum in B meson decay 



by accounting for Fermi motion of the b quark. We make the reasonable assumption 
that B meson bound state effects can be factorized from the hard subprocess as well as 
from cc fragmentation. The Fermi motion effect can be described rigorously in heavy 
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quark effective theory JTJJ, but we contend ourselves with the earlier ACCMM model 



nj. The ACCMM model is in fact consistent with the heavy quark expansion, if a 
particular relation between the b quark mass and the ACCMM model parameter p R 
is adopted (The ACCMM model then assumes a particular value for the kinetic 

energy matrix element of heavy quark effective theory.) The ACCMM model provides 
a phenomenologically viable description of energy spectra in other B decays, e.g. B — * 
Xtv t or B -> X sl . 

The basic idea of this model is quite intuitive: one imagines the b quark moving 
inside the B meson at rest with a momentum p according to some distribution with a 
width of a few hundred MeV. The cloud of gluons and light quarks in the B meson of 
the mass Mb is treated as spectator quark with mass m sp . To keep the kinematics of 
this "decay in flight" exact one introduces a so-called floating b quark mass 



m 2 b (p) = M% + ml - 2M B Jm% + p 2 . (3.21) 



The b quark is on-shell with energy E b (p) = (m 2 (p) + p 2 ) 1 ^ 2 . The b quark momentum 
distribution must be chosen ad hoc. Usually one takes a properly normalized Gaussian 
form 

$acm(p) = exp (V/pD, (3.22) 

where dpp 2 $acm(p) = 1- Implementing the kinematics of decay in flight, the J /if) 
energy distribution in the B meson rest frame (quantities without "hat" ) is then obtained 
from the spectrum in b quark decay fl3.12| ) by the convolution 

dT f , 9 ^ , , ml(p) R f dE R df 

— = / dpp 2 $ ACU (p)-^- / R . 3.23 
dE R J 2pE b (p) _ J E R dE R 

max{0,p-} E™'(p) 

The integration over the J / if) energy E R in the b quark rest frame is limited by 
pmax ■ j E R E b (p) + \p R \p m 2 (p) + M 2 R \ 

Er - mm i — — ' 2Mp) r (3 - 24) 

ET = ErE ^-} P ^. (3.25) 

The requirement E R m < E R ax leads to the following bounds on p: 

|PB ± (Mb - -Efl)] 2 - m% 
P± = 2\p R ±(M B -E R) ] ■ (3 - 26) 



The dependence of the energy spectrum (|3.23| ) on the two parameters of the ACCMM 
model, m sp and pp, is quite different. Changing the value of the spectator mass does 
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not affect the spectrum noticeably. Therefore m sp usually is fixed to 150 MeV in all 
ACCMM analyses. On the other hand the width pp of the momentum distribution 
must be chosen carefully, because the shape of the spectrum is strongly sensitive to this 
parameter. Successful fits to the lepton energy spectrum in semi-leptonic decay typically 
find p F « (300 - 450) MeV |BT|. 



D. Final result and comparison with CLEO data 

Eq. ( |3.23| ) yields the J/if> energy spectrum for B mesons decaying at rest. To compare 
with CLEO data ||16|| , we have to translate the energy spectrum ( |3.23j ) into a momentum 
spectrum 

dF E " dT (3.27) 



dp R p R dE R 

and account for the fact that B mesons have momentum ps = {Myus)/^ — M^) 1 ^ 2 »s 
482 MeV in the CLEO rest frame in which the data in |16| is presented-!^] The final 
boost from the B meson to the T(4S) rest frame is effected by 



dT p R M B f dp R dT 



dp R E R 2p B 1 Pr dp R ' 



(3.28) 



where the bounds on the J/i/j momentum 



min /n EbPr- PbE r \ , q 
p R = max <^ 0, — | , (3.29) 

. \ X^(M B ,Ml,m%) E b Pr + PbEr \ 
Pr =mm \ 2Mb ' W B ) (3 - 30) 

stem from kinematical restrictions set by the masses in the Kallen function X(x, y, z) = 
x 2 + y 2 + z 2 — 2xy — 2xz — 2yz and from the integration over the angle between the B 
and the J/ip momentum. 

Owing to the difficulties of normalizing the partial production rates discussed above 
we forsake the idea of predicting the absolute J/ip branching fraction in B decay and 
concentrate on the shape of the spectrum. We fix the absolute normalization by adjusting 
the sum of all contributions to data. This is actually equivalent to re-fitting the NRQCD 
matrix elements to data after including large higher order corrections in the velocity 



^Quantities with "tildes" refer to the CLEO or T(45) rest frame. 
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expansion. However, we do not give the result of the re-fitting, because we believe it is 
of little interest for comparison with other J/ ip production processes. 

The shape function ansatz (|2.29|) is slightly different for the different production 
channels because of the different choice of parameters (|2.30|) , (|2.31| ). Therefore the 



shape of the momentum spectrum depends somewhat on the relative contribution of the 
various channels even after adjusting the overall normalization to data. We determine 
the relative normalization of the various channels by comparing the existing information 
on the NRQCD matrix elements obtained by standard leading order NRQCD analyses. 
The colour singlet matrix element can be computed from the wave function at the origin. 



The Buchmuller-Tye potential is often adopted with the result [32 



(O^CSt)) = 9|i 5 0)l = 1.16 GeV 3 . (3.31) 

Due to our particular treatment of the colour singlet contribution as described below, 
we do not need this matrix element in B decay. The colour octet matrix elements are 
determined by fits to J/ip production in a variety of production processes. (0^( 3 Si)) 
is best determined from J/ip production in hadron-hadron collisions at large transverse 
momentum P, |33| , |12"1| , or, perhaps, from charmonium production in Z° decays |34]]. Given 



uncertainties from unknown higher order perturbative corrections a reasonable range is 

(0^( 3 Si)) = (°- 5 " L °) • 10 ~ 2 Gey3 - ( 3 - 32 ) 

The determination of the other two matrix elements from hadron-hadron collisions is 
much more uncertain. Assuming the above range for (Og^( 3 Si)) a significant constraint 
on 

M^(^ 8) , 3 P (8) ) = (O^CSo)) + — 2 (O J ^{ 3 P )) . (3.33) 

with k — 3.1 arises from the integrated J/ip branching in B decay itself 0. A reasonable 
range is 

M^CS^fP^) = (1.0 - 2.0) ■ 10- 2 GeV 3 . (3.34) 

As default we take the value M^{^ = 1.5 - 10~ 2 GeV 3 for m c = 1.5 GeV and assume that it 
originates from both parts equally. We then investigate the modification of the spectrum 
when ^Sq 8 ^, 3 Pq 8 ^) is saturated by only one of the two matrix elements and when 
the relative contribution of {Og^ ( 3 S\)) and ^S^^Pq 8 ^) is varied as allowed by the 
ranges of values given. At the end we discard the absolute normalization that would be 
implied by these values and re-fit it to data as already mentioned. 

A final comment concerns the treatment of the two-body modes B — > J/ipK and 
B — > J/ipK*, which appear as sharp resonances in the J/ip momentum spectrum. Nei- 
ther the ACCMM model nor the shape function for cc fragmentation applies to these 
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j F A = MeV: •••• 




JAF momentum (GeV/c) 



FIG. 5. Sum of colour octet modes dBR [s] /dp R [n} with n = { l sf\ 3 P {8) , 3 s{ 8) } to the 
differential branching ratio dSR,/dpn of the decay B — > J/ipX for various values of the shape 
function parameter A. The ACCMM model parameters are fixed at pp = 300 GeV and 
m s „ = 150 GeV. 



resonance contributions. Fortunately, the information provided in [16] allows us to sub- 
tract these contributions from the momentum spectrum. We then assume that the two 
resonant contributions are dual to the colour singlet contribution, while the rest of the 
spectrum corresponds to the colour octet contribution. This appears plausible, because 
we expect colour octet cc pairs to fragment into multi-body final states, with only a small 
probability that the emitted soft gluons reassemble with the spectator quark to form a 
single hadron. Hence, the experimental spectrum shown in the following plots refers to 
the CLEO data with B — * J/ipK and B — ► J/ipK* subtracted and it is compared with 
colour octet contributions only. The integrated branching fraction from the resonance 
subtracted spectrum is 0.53 %. Of course, indirect contributions from B — > ip'X and 
B — > \ C X with subsequent decay into J/ip are also subtracted. 

We have implemented the five-fold integration that leads to the final J /if) momen- 
tum spectrum into a Monte Carlo program that uses the VEGAS routine described 
in 351. Parameters are chosen as follows: G F = 1.166 • lO^GeV" 2 , |V^,| = 0.039, 
M T = 10.580 GeV, M B = 5.279 GeV and = 3.097 GeV. The Wilson coefficient 
C\g\(p) is taken at the scale fi = 4.8 GeV, which yields C[ 8 ] = 2.19. The result compared 
to data is shown in Fig. |5| for various values of the shape function parameter A (see 



the ansatz ( 2.29 ) and ( |2.31| )). Here we have fixed the ACCMM model parameters to 
Pf = 300 MeV, motivated by the CLEO analysis of semi-leptonic B decay ||31|| , and 
m sp = 150 MeV. It is clearly seen that the effect of cc fragmentation is necessary to 
reproduce the data for this choice of ACCMM parameters. Increasing A shifts the maxi- 
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0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 
JAF momentum (GeV/c) 



FIG. 6. Contributions of the different colour octet modes n = ^S^, 3 ij , 3 5"j 8 ^} to the 
sum dBR/dpR of the differential branching ratio. The shape function and the ACCMM pa- 
rameters are fixed to A = 300 MeV, p F = 300 MeV and m sp = 150 MeV. 

mum of the spectrum to lower values of pr. We get the best fit for A = 300 MeV, where 
X 2 = 30.2/20 d.o.f.. 

In order to estimate the uncertainty of this fit we investigated the sensitivity of the 
best-fit A to the variation of the relative normalization of the various cc production 
channels as described above and to the ACCMM parameter pp. Fig. ^| shows the best-fit 
result of Fig. [| broken down into the separate contributions of the three colour octet 
channels. Each channel peaks approximately at the same value pr and has similar shapes, 
although the 3 <Si contribution is somewhat broader due to the choice of c = 1.5 in ( |2.31| ). 
(Varying c between 1 and 2 does not change our fit significantly.) Thus the result of 
fitting A is rather stable under changing the weightings of the different modes. Both, 
increasing the relative contribution of and saturating it by only one of its matrix 

elements leads to variations of A of about 50 MeV. There is an obvious anti-correlation 
between the size of A and of pp, although the effect is not as large as one may expect. 
Figure [j] shows the spectra for different values of pp while A is fixed to 300 MeV. We 
obtain that the spectrum is slightly wider for higher values of pp. But even for pp = 500 
MeV the best-fit A would remain of order 200 MeV. 

We conclude from this analysis that the kinematics of soft gluon emission has to 
be accounted for to describe the data on J/ip momentum spectra and that our shape 
function model provides a satisfactory description of the spectrum shape, if the parameter 
A is chosen in the range 

A = 300:% MeV. (3.35) 
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FIG. 7. Sum of colour octet modes dBR [8] /dp R [n] with n = {^ 8) , 3 P {8) , 3 5j 8) } to the 
differential branching ratio dBR/dpn of the decay B —* J/if>X for various values of the ACCMM 
parameter pp\ The shape function parameter A = 300 MeV and the spectator mass m sp = 150 
of the ACCMM model are kept fixed. 



This result agrees perfectly with the velocity scaling rules, which lead to the estimate 
A ~ m c v 2 ~ Aq CD . It is also worth noting that the partonic spectrum behind Fig. ^| is 
a pure delta-function so that the smearing due to cc fragmentation and Fermi motion 
extends almost over the entire accessible momentum range. Only for rather small J /if) 
momentum, there would be a visible tail due to perturbative hard gluon radiation [[7] . 

Finally let us comment on the J/if> momentum spectra in based on the effect 

of Fermi motion only. (Earlier results |36| were based on the colour singlet model and 
will not be discussed.) These works also report acceptable fits of the J/if> momentum 
spectrum, however with a larger value of pp ~ 550 MeV, as one may expect when cc 
fragmentation effects are neglected. However, even this large value of pp is obtained only, 
because the K and K* resonances, which sit at large values of pn have been included, 
even though the ACCMM model cannot be applied to them. If these contributions are 
subtracted, as done in the present analysis, a satisfactory fit is not obtained with the 
ACCMM model alone. 



IV. INELASTIC J/ip PHOTOPRODUCTION 

In this section we discuss the energy spectrum in inelastic J / if> photoproduction. This is 
perhaps the most interesting application of the shape function model developed in this 
paper. The colour octet contributions to the energy spectrum have been predicted to 
increase rapidly in the endpoint region, where the J / if> approaches its maximal kinemat- 
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ically allowed energy [§,0]. If the colour octet matrix elements take the values required 
to fit the normalization of production cross sections in hadron-hadron collisions and in 
B decay, this prediction contradicts the data collected at the HERA collider ||10|| , which 



show a rather flat energy distribution. The measured distribution can in principle be de- 
scribed by colour singlet contributions alone, both at leading order and at next-to-leading 
order [^SJ in a s . 



Several solutions have been proposed to solve this problem for the NRQCD approach 
to charmonium production: 



(a) The relevant colour octet matrix elements are smaller than believed [39]. The 
colour octet contributions are always small and the shape of the energy spectrum is 
determined by the colour singlet term. 

(b) The partonic cross section is modified by intrinsic transverse momentum effects. 
Within a particular model for these effects |5J obtains a reduction of the colour octet 



cross section, while the energy dependence is essentially unmodified. 

(c) The NRQCD calculation is unreliable for large J/ip energies because of a break- 
down of the non-relativistic expansion [jnj. Resummation of the expansion as discussed 
earlier leads to folding the partonic cross section with a shape function. It is expected 
that this leads to a depression of the spectrum at large J/ip energies, because some 
energy is lost for radiation in the fragmentation of the colour octet cc pair. 

In this section we pursue suggestion (c), which has not been implemented in practice 
yet. Let us note that, irrespective of the issue of normalization, this is the only solution 
that addresses the fact that the shape of the colour octet spectrum obtained from a 
partonic calculation is unphysical for large J/ip energies. 

The section is organized as follows. In parallel with the discussion of the B decay we 
begin with kinematics and by listing the relevant partonic subprocesses 7+g — > cc[n]+g.f\ 
We then incorporate the fragmentation of the cc pair via our shape function ansatz 
and discuss the modification of the energy spectrum. For the sake of demonstration, 
we compare the result to HERA data, although we shall see that this comparison is 
problematic from a theoretical point of view. 

A. Kinematics of photoproduction 

The quantity of interest is da/dz, where 

z = (4.1) 



s Photon-quark scattering is a small correction on the scale of effects we are going to discuss, 
and relative to photon- gluon fusion. We omit these subprocesses for simplicity. 
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and p R , p p and p 7 denote the J/ip, proton and photon momentum, respectively. In the 

proton rest frame z is the fraction of the photon energy transferred to the J/ip. In the 
photon-proton centre-of-mass system (cms) we define 

P 7 = ^ (1,0, 0,-1), (4.2) 



,^(1 
' 2 v 

p R = (E R ,p T ,0,p R ), (4.4) 



^ = ^^-(1,0,0,+!), (4.3) 



where s = (p p + p 7 ) 2 is the cms energy and x g the gluon momentum fraction of the 
proton momentum. Note that z and pt refer to the physical J /if) particle. In the present 
context they cannot be identified with the corresponding quantities of the progenitor cc 
pair, which we denote by z c5 and pr,cc- Using p\ = M R , we express the J/ip energy E R 
and its longitudinal momentum p z R in terms of its transverse momentum p? and z: 

H 2z^s ' FR 2 Zy fs V ; 

The convolution with the shape function, ( |2.21| ), requires a and /3, defined by ( [2.20] ) 
in the quarkonium rest frame.** According to our convention, the i-axis is defined in 
the direction of — Pj n with P in = p 7 + p g . Writing 

Vl = (K p p ± , 0, p z 7 ) , (4.6) 

p g = (E g , -p ± , 0, f g ) , (4.7) 

P m = (E in , 0, 0, P* n ) , (4.8) 
and performing the Lorentz transformation explicitly, we obtain 

E = M n + Pr (49) 
7 2M R z ' 1 ' 

- PTZX g S 

p± \^(Mi-p* T ,x g s Z iy l4 - iUj 

z i Xg s{p 2 T -Ml) + {pl + MlY 
P ^ 2zM R X 1 l 2 {M 2 Rl -p 2 Tl x g sz 2 ) ' 1 ' 1 



** Contrary to the previous section we now use "hats" to denote quantities defined in the J/ip 
rest frame. Non-invariant quantities without hat refer to the photon-proton cms frame with z 
axis in the direction of the proton momentum. 



32 



- % < 4 - 12 > 

^ " 2M R \V*(M 2 R ,-p 2 T ,x g sz*) { • J 

with A(x, y, z) = x 2 + y 2 + z 2 — 2a;?/ — 2xz — 2yz, and 

. M| + p^ + a: g ^ 2 ~ \ l / 2 {M 2 R ^p 2 T ,x g sz 2 ) 

2M R z ' m 2M R z ' 1 j 
The previous line gives a and defined as 

a = 4 + 4-M a (3 = E in -P* n -M R (4.15) 

for given 2 and p-r of the J/^ in the cms frame. 

The Mandelstam variables that appear in the hard production amplitude for 7 + g — ► 
cc[n] + (7 are defined as 



s 



(Pg+Py) =X g S, 



t = {p c - c -p 1 )\ (4.16) 

U = (Pcc-Pg) 2 - 

We have to express them in terms of z, pt, x g and the energy ko and invariant mass k 2 
of the radiated soft partons in the J/tp rest frame. Recall that p c5 = P + l = Pr + k with 
P = (2m c ,0). Hence 

u = M cE (k) 2 -s-t, (4.17) 

where M c5 (k) 2 = M R + 2M R ko + k 2 as usual. Next parametrize the momentum of the 
cc pair by 



Pc 



(Ecc, l± cos 0, i± sin 0, 4) . (4.18) 



This introduces azimuthal angular dependence into the partonic matrix element. This 
dependence is formally small. All dependent terms are proportional to l±, and, as 
discussed in Sect. [II B 2| , such transverse momentum dependence can be neglected in the 
strict shape function limit. In our ansatz, which models the entire spectrum, we also have 
to keep the exact kinematic relations and therefore a non-trivial azimuthal average of 
the hard production amplitude appears in this case. With the help of on-shell conditions 
for the hard emitted gluon we can express the components of p c5 by 

E c - C = M R + k , 

[k 2 - a(2k - a)] 1 / 2 [/3(2A;o - 0) - k 2 } 1 ' 2 
j3 — a 

k 2 + a(3 - k (a + (3) 



1 1 = k 1 
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This result, together with the result for p 7 and a, /3, allows us to express t in terms of 
z, pt, ko, k and x g . 

Let us now turn to the hard amplitudes squared of the partonic subprocesses. We 
restrict ourselves to photon-gluon fusion, 7 + g —>■ cc[n] + g, where n represents either 
the dominant colour singlet state 3 Si or one of the colour octet states 1 So, 3 Pj, 3 Si. In 
terms of Mandelstam variables, the spin averaged expressions are |],[37],[|l|] : 

16e 2 c e 2 g 2 (2m c ) \s 2 {t + u) 2 + i 2 (u + s) 2 + u 2 (s + t) 2] 



H{ 3 S[ 1] ](s,i,u, 2m c 

H[ l sj, 8) }(s,t,u, 2m c 
H[ 3 Pf ] }(s,i,u, 2m c 



3e 2 c e 2 g 2 s su 



27 (s + t) 2 {t + u) 2 {u + s) 2 



[s + t + u) 4 + s 4 + t 4 + u 4 



(2m c )t(s + t) 2 (t + u) 2 (u + s) 



(4.20) 
(4.21) 



6e 2 c e 2 g 2 



{2m c )t{s + t) 3 {t + u) 3 (u + s) 3 



x 



F(2s 3 + 13s 2 w + 13su 2 + 2u 3 ) 
+ t 5 (4s 4 + A7s 3 u + 70s 2 u 2 + 47 su 3 + 4m 4 ) 
+ t 4 (2s 5 + 63s 4 u + 136s 3 ii 2 + 136s 2 u 3 + 63su 4 + 2u 5 ) 
+ si 3 u(A7s 4 + 132s 3 u + 190s 2 u 2 + 132su 3 + 47u 4 ) 
+ si 2 u(25s 5 + 88s 4 u + 156s 3 w 2 + 156s 2 w 3 + 88su 4 + 25m 5 ) 
+ siu(7s 6 + 38s 5 u + 78s 4 u 2 + 98s 3 u 3 + 78s 2 u 4 + 38su 5 + 7m 6 ) 



+ 7s 2 u 2 (s + u)(s 2 + su + u 2 ) 2 
^-H[ 3 S{\s,i,u,2m c ). 



H[ 3 S[ 8) ](s,t,u, 2m c 



(4.22) 
(4.23) 



Here e is the electromagnetic coupling, g s the strong coupling and e c = 2/3 the electric 
charge of the charm quark. Note that the NRQCD elements are not part of the hard 
cross sections, but included in the normalization of the radiation function $ n (/c), see 
( |2.28 ). In case of the P wave contribution, the normalization refers to (Og( 3 Po)) /m 2 
and the corresponding factor 1/m 2 is also extracted from H[ 3 Pj S ^](s, i, u, 2m c ). 

The hard amplitudes squared are then expressed as functions of z, pt, x g , k , k 2 and 
<f) and the average over the azimuthal angle <p according to ( |2.19| ) is performed. The 
double differential cross section for 'j + g^J/tp + Xis then given by 



d 2 a. 
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dp\dz 167r 2 z 



E 



a8 ?2 (/3 2 +fc 2 )/(2/3) 

dk z 



2tt 



dkn 



(a 2 +k 2 )/(2a) 



2s En{z ' P ^ X °> * ' k2) AV2(M 2 , -p 2 T , x g sz 2 ) 



(4.24) 
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The sum runs over the four cc states listed above. Note, however, that no shape function 
is required for the colour singlet contribution, since the dominant contribution to the 
colour singlet matrix element does not require emission of soft gluons. For the colour 
singlet contribution we therefore use the ordinary differential cross section on the parton 
level. The final result is obtained by folding in the gluon distribution in the proton, 
g(x g ,fj,p), and integrating over transverse momentum: 

PT.max 1 „ 

da 1P f , f , , s da, 



, J dpi JdxM*»i»)jsZ- (4 - 25) 



PT.min 



The lower integration limit for p\ usually is set by an experimental cut. In the present 
framework such a cut is needed to eliminate the contribution from the 2 — > 1 process 
7 + g — > cc[n] , smeared out over a finite range in and z by soft gluon emission in the 
fragmentation of the cc pair, and also from the initial state. The other bounds are given 
by 

p 2 T , max = (l-z)(sz-M R ), (4.26) 

x ■ = (4 27 ) 

sz(l — z) 

The minimum px cut implies that z < 1 — Pr^ m \ n / s + . . .. For large cms energy, as at the 
HERA collider, this is not a severe restriction on the z spectrum. 



B. Discussion of the energy spectrum 

The following results for the energy spectrum are obtained with the GRV94 LO gluon 
density J42| and factorization scale \if = Mr, where Mr is the J/ij) mass. We also 



use Aqqq = 0.2 GeV (consistent with GRV) and a s {Mn) = 0.275. The cms energy is 
fixed to an average photon-proton cms energy at HERA, ^Js = 100 GeV. We also choose 
m c = 1.5 GeV for the colour-singlet process. 

In Fig. | we display the J/tp energy spectrum da/dz with the J/ip transverse mo- 
mentum larger than 5 GeV. This cut is larger than the one currently used by the HERA 
experiments. However, it allows us to discuss the effect of cc pair fragmentation in a 
situation that is theoretically under better control. The curves in the upper plot of Fig. |8] 
show, as expected, that the spectrum turns over and approaches zero as z — > 1, once 
some fraction of the photon energy is lost to radiation in the fragmentation of the cc 
pair. This turn-over occurs at smaller z for larger values of the parameter A, which is 
related to the typical energy lost to radiation in the J/ip rest frame. For J/i/j production 
in B decay, we found that A « 300 MeV fitted the spectrum well. Assuming universality 
of the shape function, this is our preferred choice for photoproduction. For comparison, 
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we also display the result with A = 500 MeV. Note that these numbers refer to the J / ip 
rest frame. In another frame, such as the laboratory frame, the energy lost to "soft" 
radiation may be large, of order ErA/Mr, where Er is the J /if} energy in that frame. 

The overall normalization in Fig. |8] and the subsequent figure requires comment. 
The NRQCD matrix elements are chosen as in Sect. [Tl fJJ on B decay. As in that case 
the normalization has then to be re- adjusted to account for the fact that the effective 
charm quark mass in the hard scattering amplitude is much larger than m c = 1.5 GeV, 
conventionally assumed in fits of NRQCD matrix elements. We proceed as follows: 
The curves labelled "partonic" (total and colour singlet alone) use m c = 1.5 GeV to 
allow comparison with earlier results. For given A, and for each colour octet channel 
separately, we determine m'f defined in (|3.20| ). We then recalculate the partonic curve 



with m c = and determine a normalization ratio by dividing the result for 1.5 GeV by 
the second one in the region of z ~ 0.1 — 0.4. Finally, we compute the curve including the 
shape function with the given value of A, multiply it by this ratio and compare it to the 
partonic curve for the conventional choice m c = 1.5 GeV. The low z region is chosen to 
compute the normalization ratio, since the shape function should have little effect on the 
spectrum far away from the endpoint. As a consequence of this procedure the partonic 
result and the results for non-zero A nearly coincide for small z. The normalization 
adjustment is quite large, which reflects the strong m c dependence of the partonic cross 
sections. 

Closer inspection of the upper panel of Fig. |] shows that the spectrum for non-zero 
A increases faster for moderate z than the partonic spectrum. To understand this effect, 
we reconsider the hard amplitudes squared for the production of a colour octet cc pair in 
a 1 So or a 3 Pj state as functions of z c5 and pt,cc- For any fixed pr,cc the hard amplitudes 
squared increase as 1 / ( 1 — z cS ) 2 as2 cg — > 1 , as follows from s = —t/(l — z cS ) and u ~ — s as 
z C c — > 1. This causes the troubling increase of colour-octet contributions in the partonic 
calculation. Now, for any given z, 

z c - c = > z (4.28) 

as can be seen by going to the proton rest frame. Hence, for fixed z, the hard amplitude 
squared is evaluated at larger z c5 , when A is non-zero compared to the partonic result 
for which z c5 = z. Due to the above-mentioned behaviour of the amplitude, sampling 
the hard cross section at larger z c5 increases the spectrum. Likewise, the transverse 
momentum of the cc pair with respect to the beam axis 

Pt, c5 = (1 " zc)(x g sz c - c - M 2 C5 ) (4.29) 

differs from p\. This happens for two reasons: first, the loss of energy to radiation also 
implies a loss of transverse momentum with respect to the beam axis, if the J/if> is not 
parallel to the beam axis. Second, the J ftp can gain transverse momentum by recoil 
against the soft gluons radiated during the fragmentation process. For fixed px, this is 
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FIG. 8. The J/tp energy spectrum for \/s = 100 GeV and with a transverse momentum 
cut PT,mm = 5 GeV. Upper panel: spectrum for three values of the shape function parameter 
A = ("total partonic"), 300, 500 MeV. Dotted: colour singlet contribution alone. Lower panel: 
as upper panel but with the "modified matrix element" discussed in the text. 
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preferred to losing transverse momentum, because the production amplitude for the cc 
pair increases with smaller pt,cc- The dominant effect is the one related to z c5 > z. The 
corresponding increase of the spectrum (for A non-zero and moderate z) relative to the 
partonic spectrum is stronger as pr, mm is chosen smaller, since the hard cross section 
rises faster for smaller pr,mm (and would approach the collinear and soft divergence at 
z — 1, if PT.min = 0). Finally, at very large z, the suppression due to the radiation 
function $ n (fc) wins and turns the spectrum over to zero. 

To illustrate these remarks we define an ad hoc modification of the hard cross sections 

Cc) • 



Hr d (Zcc,PT,c 



' H n (0.9,p T ,cc) if z c5 > 0.9, p T ,cc > 1 GeV 

H n (z cc -, 1 GeV) if z C e < 0.9, p T ,cc < 1 GeV 

H n (0.9, 1 GeV) if z c - c > 0.9, p T , C c < 1 GeV 

H n (zec,PT,cc) else 



(4.30) 



The energy spectra analogous to the upper panel of Fig. |] but with hard cross sections 
modified in this way are shown in the lower panel of this figure. The partonic cross 
section is modified only for z > 0.9 by construction. The spectra for non-zero A are 
reduced already at smaller z, which shows the sensitivity to z cg > 0.9 at such small z. 

We emphasize that no physical significance should be attached to the lower panel of 
Fig. |8|. The growth of the colour octet cross sections at large z is physical and reflects 
the growth of 2 — > 2 cross sections at large rapidity difference due to t-channel gluon 
exchange. In the endpoint region t ~ — Pxmm an d s ~ — u ~ Vrmm/^~ z ) so that s \t\. 
Higher order corrections to the spectrum would involve logarithms of s/ (— t). Summation 
of these logarithms with the BFKL (Balitsky-Fadin-Kuraev-Lipatov) equation increases 
the parton cross section in the endpoint region. 

After this discussion for large transverse momentum of the J/if>, we display the result 
for the energy spectrum with the additional requirement > 1 GeV, which we compare 
to data from the HI and ZEUS collaborations [Kj. Fig. ^ shows again the conventional 
partonic calculation compared to the calculation with two values of A. The lower panel 
refers to the ad hoc modification of the hard cross sections according to ( J4.3U ) . 



The qualitative features evident in the upper panel follow from the previous discus- 
sion. At large z the spectrum turns over, but at intermediate z, including the entire 
region where data exist, there is a large enhancement, which follows from the fact that 
the partonic matrix element is sampled very close to z C c = 1. Taken at face value, the 
disagreement with data is worse after accounting for cc fragmentation effects. However, 
the theoretical prediction with small transverse momentum cut is unreliable at large z. 
With no pt cut at all, we expect that the z spectrum is drastically modified at large 
z after accounting for the 2 —>■ 1 process, the virtual corrections to it, and soft gluon 
radiation from the initial gluon. Owing to the sensitivity to large z C c, the theoretical pre- 
diction is more sensitive to these modifications when gluon radiation in cc fragmentation 
is included. An indication of this is provided by plotting the spectrum with the modified 
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FIG. 9. The J ftp energy spectrum at y/s = 100 GeV and with pr > 1 GeV compared to 
HERA data fTO] . Upper and lower panel as in Fig. ||. Solid (dash-dotted, dashed) lines refer 
to A = 300 (500, 0) MeV. Dotted: colour singlet spectrum alone. 
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partonic cross section. This modification of the partonic cross section, although ad hoc, 
may give a qualitative representation of the effects to be expected from soft gluon re- 
summation. The lower panel of Fig. ^ shows that the unphysical enhancement is largely 
reduced in this case, although it does not disappear completely. If reality turned out to 
resemble the lower panel, it would be difficult to disentangle colour octet contributions, 
given the additional normalization uncertainties of both, the colour singlet and the colour 
octet contributions. In this case a J/f polarization measurement would provide useful 
additional information |37| . 



The results of this analysis can be summarized as follows: with the small transverse 
momentum cut on the J/ip currently used by both HERA collaborations, the region 
z > 0.7 is beyond theoretical control. This remains true even after resummation of large 
higher order NRQCD corrections via the shape function, since the hard partonic cross 
section is sensitive to other modifications that are also difficult to control theoretically 
at such small transverse momentum. At present, the experimental data cannot be in- 
terpreted as providing evidence for or against the presence of colour octet contributions 
in photoproduction. It is not necessary to reduce the colour octet matrix elements as 
suggested in to arrive at this conclusion. This is welcome as matrix elements of 
the order quoted in Sect. fT| seem to be needed to account for the observed branching 
fraction of B —>■ J/ipX. 

The situation in photoproduction remains unsatisfactory. In our opinion, nothing is 
learnt on quarkonium production mechanisms, if a small transverse momentum cut is 
used. We therefore recommend that future increases in integrated luminosity should not 
be used to reduce the experimental errors on the present analysis, but to increase the 
transverse momentum cut at the expense of statistics. 



V. CONCLUSION 

In this paper we provided a first investigation of the kinematic effect of soft emission in 
the fragmentation process cc[n] — > J/ip + X. In the NRQCD factorization approach to 
inclusive quarkonium production these effects appear as kinematically enhanced higher 
order corrections in the NRQCD expansion [|l7],[il|, which become important near the 
upper endpoint of quarkonium energy /momentum spectra. The shape function formal- 
ism discussed in fT7| , [l"8| resums these corrections and allows us to extend to validity of the 
NRQCD approach closer to the endpoint, although the entire spectrum is not covered 
even after this resummation. In the present paper, we implemented the kinematics of soft 
gluon radiation exactly and used an ansatz for the probability of radiation of soft gluons. 
This allows us to cover the entire energy spectrum, although in a model-dependent way. 
The model is consistent with the NRQCD shape function formalism in the energy region 
where the later applies. This situation is similar to the relation of the ACCMM model 
to the heavy quark expansion in inclusive semi-leptonic decays of B mesons. The main 
result is provided by ( 2.21| ), which applies to a general inclusive quarkonium production 
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process, when the partonic final state before fragmentation consists of a cc pair and one 
additional massless, energetic parton. 

We then proceeded to two applications of the formalism. These applications are not 
necessarily the simplest ones conceivable, but they seem to be most interesting. We first 
considered J /if) production in B decay, which proceeds through colour octet states by 
a large fraction. In this case the effect of emission in fragmentation of colour octet cc 
pairs has to be disentangled from Fermi motion of the b quark in the B meson. We 
found that the description of the spectrum improves significantly, when soft radiation is 
included, and if the parameter A that controls the energy scale for soft radiation is chosen 
around 300 MeV. The shape function defined in |L8j is production process independent. 
Assuming universality of our ansatz over the entire energy range, the same ansatz is 
used for inelastic J / ip photoproduction. We found that the energy spectrum turns over 
at z pa 0.8 — 0.9, to be compared with the partonic spectrum that rises towards z — 1. 
However, at z < 0.8, the colour octet contributions to the spectrum still grow faster than 
the colour singlet contribution. Due to the increase of the partonic cross section, the 
increase is in fact faster in this intermediate z region after cc fragmentation effects are 
included. We also concluded that the transverse momentum cut pr > 1 GeV, presently 
used by the HERA experiments, is too small to arrive at a reliable theoretical prediction. 
Hence, no conclusion regarding colour octet contributions and the validity of the NRQCD 
formalism can presently be drawn from HERA data. 

The formalism developed in this paper could be applied to other production processes, 
in which the J/ip energy is measured. Another interesting extension is quarkonium 
decays, when the energy of one of the decay particles is measured, such as the photon 
energy in r\ c — > 7 + X and J/ip — > 7 + X. Since decay processes are less affected by 
theoretical uncertainties related to colour recombination and initial state radiation than 
production processes, this may lead to theoretically better controlled applications of the 
shape function formalism. 
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